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Abstract 


Smoothed  particle  hydrodynamics  (SPH)  is  a  method  of  continuum 
mechanics  analysis  in  which  materials  are  modeled  by  a  discrete  set  of 
particles.  The  SPH  code  Magi  has  been  used  to  simulate  the  penetra¬ 
tion  of  a  semi-infinite  steel  target  by  tungsten  alloy  disks  traveling  at  an 
initial  impact  velocity  of  2  km/s.  Calculations  were  performed  to  sim¬ 
ulate  experimental  configurations  using  one,  two,  and  four  disks.  All 
disks  had  a  constant  length-to-diameter  ratio  of  0.125.  The  computed 
penetration  depth  into  the  target  material  is  compared  to  the  exper¬ 
imental  data  for  each  case.  The  study  included  a  set  of  calculations 
in  which  the  problem  resolution  was  varied  to  determine  the  ability 
of  the  method  to  converge  on  a  penetration  depth  as  the  number  of 
particles  in  the  problem  was  increased.  Advantages  and  limitations  ob¬ 
served  in  the  application  of  SPH  to  the  field  of  penetration  mechanics 
are  discussed. 
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1.  INTRODUCTION 


Penetration  and  perforation  of  solids  have  long  been  of  interest  to  the  military.  The  de¬ 
velopment  of  armors  to  improve  the  survivability  of  military  vehicles,  the  design  of  advanced 
penetrators  and  warheads  to  increase  the  lethality  of  weapon  systems,  and  the  application  of 
blast  and  fragmentation  mitigation  devices  to  vehicles  and  buildings  are  fields  of  study  requir¬ 
ing  knowledge  of  penetration  mechanics.  The  safety  of  nuclear  reactor  containment  vessels, 
the  development  of  lightweight  body  armors  for  police  forces,  the  protection  of  spacecraft 
from  meteoroid  impact,  and  anti-terrorist  technology  are  but  a  few  non-military  applications 
of  the  study  of  penetration  phenomena.  A  thorough  review  of  the  fundamentals  of  penetra¬ 
tion  and  perforation  of  solids  and  their  application  to  practical  problems  has  been  prepared 
by  Goldsmith  (1960),  Johnson  (1972),  Backman  and  Goldsmith  (1978),  and  Zukas  et  al. 
(1982,  1990). 

Impact  and  penetration  events  typically  involve  the  interaction  of  vaxious  types  of  materi¬ 
als,  the  propagation  of  strong  shock  waves,  large  deformations  and  fracture  of  solid  materials, 
high  velocities,  large  strain  rates,  and  detonation  of  explosives.  Numerical  methods  in  con¬ 
tinuum  mechanics  have  been  used  to  model  penetration  phenomena  since  the  1950s.  Early 
applications  lacked  the  ability  to  model  the  strength  characteristics  of  solids  and  considered 
only  the  hydrodynamic  (i.e.,  density  dominated)  phenomena.  Material  strength  and  fracture 
criteria  were  implemented  in  these  “hydrocodes”  during  the  1970s.  However,  the  develop¬ 
ment  of  advanced  material  models  is  still  an  active  area  of  research  today,  with  particular 
emphasis  on  material  damage,  in-homogeneous  materials,  ceramics,  and  non-isotropic  mate¬ 
rials.  A  concise  review  of  material  behavior  at  high  strain  rates  has  been  edited  by  Blazynski 
(1987). 

Most  numerical  methods  for  continuum  mechanics  analysis  involve  the  use  of  a  grid  or 
mesh  to  discretize  the  computational  field  (also  referred  to  as  the  computational  domain). 
For  each  cell  or  element  within  the  computational  field,  the  conservation  equations  are  ap¬ 
plied  on  a  control  volume  basis.  The  global  partial  differential  equations  representing  the 
conservation  variables  are  transformed  into  algebraic  equations  that  are  applied  locally  at 
each  cell.  Numerical  methods  of  this  type  fall  into  two  categories:  Eulerian  and  Lagrangian. 
In  the  Eulerian  formulation,  the  conservation  equations  are  applied  from  the  perspective  of 
events  taking  place  at  a  fixed  point  in  space.  This  formulation  leads  to  a  computational 
method  in  which  the  grid  is  fixed,  and  material  motion  is  accomplished  through  transport  of 
the  material  across  the  boundaries  of  cells  within  the  grid.  In  the  Lagrangian  formulation, 
the  conservation  equations  are  derived  on  the  basis  of  an  observer  that  moves  with  the  mate¬ 
rial.  Consequently,  in  Lagrangian  methods,  the  grid  moves  with  the  material  and  there  is  no 
transport  of  material  across  cell  or  element  boundaries.  Previous  reviews  of  one-dimensional 
(ID)  and  two-dimensional  (2D)  codes  for  wave  propagation  and  impact  have  been  prepared 
by  Mescall  (1974),  Von  Rieseman  et  al.  (1974),  Herrman  (1975),  and  Belytschko  (1975).  A 
review  of  three-dimensional  (3D)  wave  propagation  codes  for  impact  problems  is  given  by 
Zukas  et  al.  (1982,  1990).  Benson  (1992)  recently  documented  a  comprehsive  review  of  the 
physics  and  numerical  methods  for  modeling  the  dynamics  of  impacting  solids. 
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Each  of  these  two  approaches  has  a  unique  set  of  strengths  and  weaknesses  when  applied 
to  the  study  of  impact,  shock  wave  propagation,  large  deformations,  and  material  fracture.  In 
Lagrangian  methods,  the  material  interfaces  are  defined  by  the  shape  of  the  mesh.  Material 
interfaces,  by  definition,  occur  naturally  along  cell  interfaces.  Because  the  mesh  moves  with 
the  material,  time-dependent  effects  on  the  material  are  properly  captured.  However,  for  3D 
cases  involving  complex  geometries,  grid  generation  can  be  a  time-consuming  task.  Where 
large  material  deformations  occur,  Lagrangian  methods  can  experience  regions  where  the 
mesh  is  highly  distorted,  which  can  result  in  numerical  instabilities  and  a  very  small  time 
step.  In  regions  where  fracture  is  involved,  it  may  be  necessary  for  Lagrangian  methods  to 
automatically  redefine  the  mesh  to  capture  dislocations  and  spallation  of  material. 

In  Eulerian  methods,  relatively  simple  structured  grids  are  typically  employed  and  the 
material  geometries  are  defined  by  the  insertion  of  material  into  the  grid.  As  a  result,  grid 
generation  for  Eulerian  methods  can  be  a  trivial  task.  Because  the  grid  is  fixed,  Eulerian 
methods  are  not  prone  to  mesh  distortions  in  cases  involving  large  deformations.  Material 
motion  by  transport  across  cell  boundaries  allows  for  fracture  and  spallation  without  the 
need  for  redefinition  of  the  grid.  The  primary  disadvantage  of  Eulerian  methods  lies  in  their 
inability  to  track  material  interfaces.  Regions  of  the  grid  where  different  materials  are  in 
contact  typically  result  in  computational  cells  containing  multiple  materials  and  possibly 
some  void.  Material  failure  is  typically  handled  by  the  insertion  of  void  into  cells  that  meet 
some  material  failure  criteria.  Computation  of  a  thermodynamic  state  and  application  of 
material  models  for  such  mixed  cells  usually  involves  iterative  computations  weighted  by  the 
relative  masses  or  volumes  of  the  materials  in  the  cell.  Material  interfaces  in  such  cases  must 
be  generated  artificially  based  on  an  interface  tracking  method.  Material  transport  across 
cell  boundaries  may  also  result  in  loss  of  time-dependent  effects  in  material  characteristics. 

These  weaknesses  of  grid-based  Euler  and  Lagrange  methods  provide  sufficient  moti¬ 
vation  to  study  the  application  of  alternate  methods  for  modeling  penetration  mechanics 
problems.  This  report  documents  the  application  of  the  smoothed  particle  hydrodynamics 
(SPH)  method  for  simulating  high  velocity  impact  experiments  performed  at  the  U.S.  Army 
Research  Laboratory  (ARL). 


2.  SMOOTHED  PARTICLE  HYDRODYNAMICS 

SPH  is  a  meshless  Lagrangian  particle  method  that  uses  interpolation  theory  to  establish 
estimates  of  field  variables  at  a  point.  The  Lagrangian  formulation  results  in  an  accurate 
representation  of  material  shapes  and  interfaces.  The  absence  of  an  underlying  mesh  elimi¬ 
nates  the  possibility  of  time  step  instabilities  that  result  from  severe  mesh  distortion.  Unlike 
Eulerian  methods  that  require  grid  in  regions  of  interest  that  may  not  contain  any  material, 
SPH  employs  particles  only  where  material  exists,  resulting  in  a  natural  treatment  of  voids. 

The  SPH  code  that  forms  the  basis  of  the  current  study.  Magi,  (Libersky  et  al.,  1993; 
Randles  &  Libersky,  1996)  employs  an  explicit  time-marching  algorithm  to  advance  the 
solution  from  time  t  to  time  t  5t.  Particle  properties  of  density,  velocity,  energy,  deviatoric 
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stress,  and  position  are  computed  for  time  t  +  5t,  based  on  values  at  time  t  according  to 
Equations  1  through  5,  respectively. 


pn+l  =  (1  _  DSe) 

(1) 

Un+1/2  ^  jjn-1/2  ^  1  ^  p 

(2) 

En+l  =  ^ 

(3) 

s:v  = 

(4) 

=xl  + 

(5) 

In  these  equations,  the  n  + 1  superscript  represents  the  solution  at  the  new  time,  t  + 
while  the  n  superscript  identifies  data  from  the  time  of  the  previous  integration  cycle,  t.  The 
time  step  identified  by  5t^  represents  the  time  step  that  is  used  to  advance  the  solution  to 
the  new  time,  while  identifies  the  time  step  used  in  the  previous  integration  cycle.  The 
time  step  is  computed  by  finding  the  minimum  u;h/(c  +  s)  over  all  particles,  in  which  c  is  the 
sound  speed,  s  is  the  particle  velocity,  h  is  the  smoothing  length  used  in  the  interpolation 
functions,  and  w  is  &  constant  stability  factor  between  0  and  1  (a  value  of  0.3  is  typical). 
The  subscripts  a  and  (3  in  Equations  1  through  5  denote  spatial  directions  and  are  used 
individually  as  vector  notation  for  positions  and  velocities  and  together  to  identify  stress 
tensors. 

The  variables  D,  F,  G,  and  H  represent  the  volumetric  strain,  total  acceleration,  work 
per  unit  mass,  and  stress  rate  on  a  particle,  respectively.  These  values  are  computed  from 
the  interpolation  function,  resulting  in  “kernel  estimates”  of  field  variables  at  a  point.  To 
compute  the  new  values  of  the  field  variables  for  a  particle,  that  particle’s  previous  values 
and  the  values  of  its  neighboring  particles  are  used  with  the  B-spline  smoothing  function 
shown  in  Equation  6: 

f  ^  (I  -  0  <  <  1 

=  <  ^  (2  - l<u  <2  (6) 

0  1/  >  2 

\ 

in  which  u  —  |x,-  —  Xj\/h.  This  smoothing  function  generates  a  weighting  factor  for  the 
values  of  neighboring  particles,  which  is  based  on  the  smoothing  length,  h,  and  the  absolute 
distance  of  the  neighboring  particle  from  the  particle  of  interest,  |Xj  — Xj|.  This  interpolation 
function  is  always  positive,  maximum  at  a  distance  of  zero,  and  reaches  zero  at  a  distance 
of  2h.  Thus,  to  compute  the  new  interpolation  functions  for  a  particle,  the  only  neighboring 
particles  that  have  an  influence  are  those  within  a  distance  of  two  smoothing  lengths  of  the 
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Figure  1.  The  interpolation  smoothing  function. 

particle  of  interest.  Figure  1  shows  a  plot  of  the  B-spline  as  a  function  of  relative  distance 
from  the  particle  of  interest. 

Because  SPH  is  a  meshless  method,  there  is  no  underlying  connectivity  logic  for  particle 
interactions  for  a  particular  computation.  The  set  of  particles  whose  values  must  be  used  to 
compute  the  new  state  variables  for  a  given  particle  depends  on  the  positions  of  the  particles. 
Because  the  particle  positions  are  updated  at  every  time  step,  a  list  of  neighboring  particles 
must  be  generated  for  each  particle  for  every  integration  cycle.  The  set  of  neighboring  par¬ 
ticles  that  may  interact  with  a  given  particle  is  defined  by  the  smoothing  function.  Because 
the  smoothing  function  depends  on  the  smoothing  length,  /i,  a  small  smoothing  length  will 
result  in  fewer  particles  in  the  neighbor  lists  (and  vice  versa). 

Thus,  for  every  iteration  in  time,  each  particle  must  compare  its  position  to  all  other 
particles  in  the  computation  and  must  build  a  neighbor  list  before  the  state  variables  can  be 
updated.  For  a  calculation  containing  N  particles,  this  implies  a  workload  of  position 
comparisons  for  every  iteration  in  time.  To  decrease  the  workload  associated  with  finding 
neighbors,  an  artificial,  rectilinear  grid  is  generated  over  the  computational  space.  Particles 
are  then  assigned  to  the  cells  of  this  artificial  grid.  The  cells  within  the  artificial  grid  have 
sides  of  length  2h,  the  range  of  the  smoothing  function.  As  previously  stated,  the  time¬ 
marching  scheme  is  explicit.  Thus,  for  a  particular  time  step,  a  particle  may  only  interact 
with,  at  most,  the  other  particles  in  its  artificial  cell  and  those  in  the  immediate  neighbor  cells. 
This  approach  reduces  the  workload  of  constructing  neighbor  lists  from  N'^  to  approximately 
Nlog{N). 
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3.  SEGMENTED  PENETRATOR  EXPERIMENTS 


The  penetration  performance  of  high  density,  low  length-to-diameter  (L/D)  ratio  pro¬ 
jectiles  impacting  steel  targets  has  been  a  topic  of  considerable  interest  in  penetration  me¬ 
chanics  because  of  the  speculated  performance  of  a  segmented  rod  projectile.  This  has  been 
spurred  by  the  observation  that  relative  projectile  performance,  i.e.,  penetration  per  unit 
length  (P/L),  increases  as  L/D  decreases,  provided  that  the  impact  velocity  is  relatively 
high.  Computational  and  experimental  research  to  date  has  focused  on  penetrators  shaped 
as  either  right-circular  cylinders  with  an  L/D  of  one  or  slightly  greater,  or  spheres  (Bris- 
senden,  1992;  Charters  et  al.,  1990;  Charters,  1987;  Charters  k  Orphal,  1990;  Cline  et  al, 
1989;  Cuadros,  1990;  Frank  k  Zook,  1990;  Hauver  k  Melani,  1990;  Hohler  k  Stilp,  1990; 
Holland  et  al.  1990;  Hunkier,  1989;  Kivity  et  al.  1989;  Naz  k  Lehr,  1990;  Raatschen  et  al., 
1989;  Scheffler  k  Zukas,  1990;  Scheffler,  1990;  Sorensen  et  al.,  1991;  Tate,  1990;  Zukas, 
1990).  De  Rosset  and  Sherrick  (1996)  modeled  segmented  rod  performance  at  ordnance 
velocity  for  high  density  tungsten  alloy  segments  with  an  L/D  of  one.  They  observed  that 
multiple-segment  rod  performance  was  less  than  that  predicted  by  simply  multiplying  single 
segment  performance  by  the  total  number  of  segments,  because  of  interactions  with  residual 
segment  material  in  the  penetration  cavity. 

Recently,  computational  and  experimental  studies  have  focused  on  characterization  and 
understanding  of  the  penetration  mechanics  of  high  density  metallic  projectiles  with  an  L/D 
of  less  than  one  (Bjerke  et  al.,  1992;  Orphal  et  al.,  1990;  Orphal  et  al.,  1992).  Herbette 
(1989)  noted  a  dramatic  increase  in  P/L  for  steel  disks  with  an  L/D  ratio  of  1/30  impacting 
aluminum  targets  at  2  km/s,  when  compared  to  penetrators  with  considerably  greater  L/D. 
Orphal  and  Fransen  (1990)  also  reported  a  significant  increase  in  P/L  as  projectile  L/D  was 
reduced  from  1  to  1/8  for  tungsten,  tungsten  alloy,  and  tantalum  alloy  projectiles  impacting 
steel  targets  at  striking  velocities  between  1.5  km/s  and  7.5  km/s. 

Experiments  have  been  conducted  at  ARL  in  which  tungsten  alloy  disk  penetrators 
were  launched  into  semi-infinite  rolled  homogeneous  armor  (RHA)  steel  targets  at  a  striking 
velocity  of  2.0  km/s  (Bjerke,  1997).  The  disks  were  1.27  cm  in  diameter  and  had  a  thickness 
of  0.3175  cm  (L/D  =  1/8).  Experiments  were  performed  using  one,  two,  and  four  disks.  For 
the  cases  using  more  than  one  disk,  the  spacing  between  adjacent  disk  faces  was  5.08  cm 
(two  disk  diameters).  The  experiments  were  designed  to  have  all  the  disks  aligned  axially 
with  the  faces  of  the  disks  parallel  to  the  face  of  the  target  material  on  impact.  A  schematic 
diagram  of  the  experimental  configuration  is  provided  in  Figure  2. 


4.  TWO-DIMENSIONAL  AXISYMMETRIC  SIMULATIONS 

A  series  of  calculations  was  performed  in  a  2D  axisymmetric  geometry  to  simulate  the 
experiments  conducted  by  Bjerke.  The  simulation  of  the  tests  in  an  axisymmetric  geometry 
is  reasonable  as  the  tests  resulted  in  nearly  perfect  axial  alignment  of  the  disks  and  normal 
impact  of  the  disks  with  the  target.  The  semi-infinite  RHA  target  material  was  modeled  as 
a  continuous  set  of  particles  with  a  radius  of  10  cm  and  a  depth  of  10  cm.  Transmissive 
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Tungsten  Disk 
D  =  2.54  cm 
L/D  =  1/8 
Vel  =  2 .0  km/s 


5.08  cm 
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boundaries  were  placed  at  the  outer  radius  of  the  target  and  at  the  bottom  of  the  target 
(opposite  the  point  of  impact)  to  model  the  semi-infinite  nature  of  the  target.  The  hydro- 
dynamic  behavior  of  the  tungsten  alloy  disks  and  the  R,HA  target  materials  was  modeled 
with  a  linear  Hugoniot  shock-particle  velocity  equation  of  state.  An  elastic/perfectly  plastic 
plasticity  model  was  used  for  all  materials.  A  summary  of  the  material  properties  is  provided 
in  Table  1. 


Table  1.  Material  Properties. 


Tungsten 

RHA 

Density  (g/cm®) 

17.270 

7.842 

Sound  Speed  (m/s) 

4237 

4529 

Slope  of  Us-Up 

1.24 

1.50 

Griineisen  Parameter 

1.80 

1.84 

Specific  Heat  at  Constant  Volume  (cal/g-K) 

0.036 

0.107 

Yield  Strength  (kbar) 

19.3 

7.0 

Shear  Modulus  (Mbar) 

1.606 

0.816 

Fracture  Pressure  (kbar) 

-30.0 

-25.0 

4.1.  Single-Disk  Simulations 

The  first  SPH  calculation  in  the  study  simulated  the  impact  of  a  single  disk.  This 
calculation  employed  a  smoothing  length,  h,  of  0.0635  cm,  with  one  particle  per  smoothing 
length.  This  resulted  in  a  distribution  of  five  particles  across  the  thickness  of  the  disk  and  20 
particles  across  the  radius,  for  a  total  of  100  particles  in  the  disk.  The  target  was  modeled 
with  the  same  smoothing  length  characteristics  and  consequently  had  identical  particle  sizes 
as  the  disk.  In  all,  the  target  contained  24,649  particles.  As  shown  in  Figure  2,  the  disk  was 
initially  positioned  above  the  target  with  its  impact  face  positioned  1  mm  from  the  impact 
face  of  the  target.  The  calculation  was  run  to  a  simulated  time  of  100  ^s,  requiring  2,527 
time  iterations. 

When  the  disk  strikes  the  target,  one  compression  wave  is  established  in  the  target  and 
another  is  established  in  the  disk.  Because  the  disk  is  very  thin  (0.3175  cm),  the  initial 
compression  wave  propagating  through  it  reaches  the  back  face  of  the  disk  at  approximately 
0.75  yus,  at  which  time,  it  reflects  back  toward  the  impact  face  of  the  disk  in  the  form  of  an 
expansion  wave.  This  expansion  wave  propagation  causes  the  disk  material  to  decelerate, 
resulting  in  the  formation  of  a  void  region  between  the  disk  and  the  target  at  the  centerline 
of  the  impact.  This  void  grows  and  then  ultimately  collapses,  resulting  in  renewed  contact 
between  the  disk  and  target  material  at  the  centerline  as  the  penetration  process  continues. 

This  process  is  illustrated  in  Figure  3  which  provides  a  set  of  material  plots  at  times  of 
4  yus,  16  yus,  27  fj,s,  and  100  yus  during  the  simulation.  Since  the  simulation  was  performed 
using  symmetry  about  the  y-axis,  the  materials  in  these  plots  have  been  reflected  about  the 
y-axis  to  show  the  actual  physical  geometry.  While  the  target  material  was  modeled  having 
a  radius  of  10  cm,  the  plots  are  limited  to  a  radius  of  3  cm.  The  plot  at  4  (is  shows  the  initial 
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penetration  of  the  disk  into  the  target  material.  The  event  at  this  time  is  characterized  by 
the  nearly  horizontal  material  interface  and  the  beginning  of  the  formation  of  a  void  region 
at  the  disk-target  interface.  At  16  yus,  this  void  region  has  grown  to  its  maximum  and  then 
begins  to  collapse  until  the  two  materials  are  again  in  contact  at  the  centerline  at  27  fis.  The 
material  plot  at  100  j^s  shows  the  final  shape  of  the  penetrator  and  the  target.  The  predicted 
depth  of  penetration  is  1.34  cm,  for  a  P/L  of  4.2.  The  predicted  P/L  is  12.4%  lower  than 
the  experimental  mean  and  is  also  in  good  agreement  with  previous  computational  values 
(Bjerke  et  ah,  1992). 

Two  Lagrangian  tracers  were  employed  in  this  calculation  to  gather  particle  history 
data.  One  tracer  was  attached  to  the  centerline  particle  on  the  impact  face  of  the  disk. 
The  other  tracer  was  placed  in  the  target  material,  attached  to  the  centerline  particle  at  the 
impact  face.  Analysis  of  the  position  and  velocity  histories  of  these  two  particles  reveals  the 
growth  and  collapse  of  the  void  region  between  the  two  materials.  Figure  4  provides  the 
y-position  (axial)  histories  of  these  two  particles  as  well  as  a  history  of  the  distance  between 
the  two  particles.  Figure  5  provides  the  y-velocity  (axial)  histories  of  the  tracers.  The  data 
provided  by  the  tracers  yield  the  center  position  of  the  particles.  Because  each  particle  has  a 
finite  radius  equivalent  to  half  the  smoothing  length,  h,  the  minimum  difference  between  two 
adjacent  particles  that  are  in  contact  is  approximately  h.  This  is  illustrated  in  the  position- 
difference  curve  in  Figure  4.  Initially,  the  disk  tracer  particle  moves  downward  and  as  the 
impact  process  begins,  the  target  tracer  particle  is  accelerated.  At  approximately  2  /^s,  the 
distance  between  the  two  tracers  is  at  a  minimum. 

After  the  initial  contact  between  the  disk  and  the  target,  propagation  of  the  expansion 
wave  through  the  disk,  coupled  with  concentration  of  radial  waves  at  the  centerline,  deceler¬ 
ates  the  disk  particles  near  the  centerline  and  causes  the  formation  of  the  void  region  between 
the  disk  and  the  target.  This  void  region  grows  to  a  maximum  at  16  fis  and  then  decreases 
until  the  two  tracers  are  again  in  contact  at  27  fiS,  as  illustrated  by  the  position-difference 
curve  in  Figure  4. 

The  velocity  histories  in  Figure  5  show  the  initial  velocity  of  the  disk  tracer  as  -2  km/s, 
the  impact  velocity  of  the  disk.  The  velocity  of  the  target  tracer  is  initially  zero  and  is 
accelerated  until  the  two  tracers  have  a  common  velocity  of  approximately  -1.25  km/s  at 
approximately  1  /us.  The  two  tracers  share  a  common  velocity  until  approximately  2.5  /us, 
after  which,  the  velocity  of  the  target  tracer  is  larger  than  that  of  the  disk  tracer,  resulting 
in  the  formation  of  the  void  region.  By  16  /j,s,  the  void  has  reached  its  maximum  and  is 
marked  by  the  crossing  of  the  two  velocity  histories.  After  this  time,  the  void  closes  until  the 
two  tracer  particles  are  again  in  contact  and  share  a  common  velocity  at  30  fis.  By  38  fis, 
the  maximum  penetration  depth  has  been  achieved. 

After  this  baseline  single-disk  simulation  was  completed,  three  additional  calculations 
were  performed  to  assess  the  ability  of  the  code  to  converge  on  a  solution  with  increasing 
resolution.  As  stated  earlier,  the  baseline  calculation  employed  a  particle  size  of  0.0635  cm 
to  produce  five  particles  across  the  thickness  of  the  disk.  The  next  three  single-disk  cases 
employed  particle  sizes  of  0.03175  cm,  0.02117  cm,  and  0.015875  cm,  which  resulted  in  10, 15, 
and  20  particles  across  the  disk  thickness,  respectively.  Each  of  these  calculations  was  run 
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Figure  3.  Single-disk  impact  at  4  fis,  16  jus,  27  /^s,  and  100  ;us  (/i=0.0635  cm). 
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Figure  4.  Tracer  position  history  for  single-disk  impact  (/i=0.0635  cm). 


Figure  5.  Tracer  velocity  history  of  single-disk  impact  (/i=0.0635  cm). 
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for  a  simulated  time  of  100  fis.  Material  plots  for  each  of  the  four  calculations,  showing  the 
final  penetration  into  the  target,  are  provided  in  Figure  6.  Each  of  these  results  shows  the 
target  to  have  a  final  crater  shape  that  is  approximately  spherical,  with  a  small  additional 
depression  near  the  centerline.  These  plots  also  show  some  erosion  of  the  penetrator  and  a 
similar  depressed  region  near  the  centerline. 

The  plot  in  Figure  7  shows  the  final  penetration  depth  for  each  single-disk  calculation 
as  a  function  of  the  particle  size.  This  figure  shows  that  the  penetration  depth  for  the 
baseline  case  is  much  less  than  that  for  the  higher  resolution  cases  and  that  the  penetration 
depth  converges  as  the  number  of  particles  is  increased.  The  case  that  employed  five  particles 
across  the  disk  thickness  produced  a  final  penetration  depth  of  1.34  cm,  while  the  cases  using 
10,  15,  and  20  particles  all  produced  penetration  depths  of  1.43  cm.  These  results  indicate 
that  the  use  of  five  particles  to  model  a  critical  dimension  (in  this  case,  the  disk  thickness)  is 
insufficient  to  properly  capture  the  physics  of  the  penetration  event.  Furthermore,  the  results 
show  that  there  is  no  benefit  in  using  more  than  10  particles  across  the  critical  dimension 
for  this  problem. 

Also  included  in  Figure  7  is  a  data  marker  to  represent  the  experimental  data  collected 
for  the  single-disk  configuration  (Bjerke  et  ah,  1992).  A  total  of  four  tests  was  performed,  re¬ 
sulting  in  a  mean  penetration  depth  of  1.53  cm,  with  a  sample  standard  deviation  of  0.26  cm. 
The  experimental  data  are  represented  in  Figure  7  as  a  data  point  at  the  mean  value  with 
an  error  bar  whose  size  is  one  standard  deviation.  The  figure  shows  that  the  computational 
result  using  five  particles  across  the  disk  thickness  falls  below  the  experimental  data  range, 
while  the  higher  resolution  calculations  fall  within  the  variation  of  the  experimental  data, 
with  a  final  penetration  depth  that  is  6.5%  less  than  the  experimental  mean. 

These  findings  are  consistent  with  the  general  application  of  Eulerian  methods  to  prob¬ 
lems  involving  strong  shock  wave  propagation.  In  these  cases,  it  is  considered  good  practice 
to  employ  at  least  five  to  ten  Eulerian  cells  across  the  smallest  critical  dimension  in  a  simu¬ 
lation.  As  in  the  use  of  Eulerian  codes,  the  implementation  of  this  practice  in  SPH  is  strictly 
dictated  by  the  overall  memory  and  processor  limitations  of  the  computer  being  used  to  run 
the  simulation.  Improved  problem  definition  will  dramatically  increase  the  memory  and  run 
time  requirements  since  the  workload  in  a  2D  simulation  is  proportional  to  while 

that  of  a  3D  simulation  is  proportional  to  {llhy.  To  illustrate  this  point,  the  number  of 
particles,  number  of  time  integration  cycles,  and  the  computer  time  are  summarized  for  each 
single-disk  calculation  in  Table  2. 


4.2.  Multiple-Disk  Simulations 

Additional  disk  impact  experiments  were  performed  in  which  sets  of  two  or  four  disks 
were  launched  into  the  semi-infinite  RHA  target  material.  For  all  of  the  multiple-disk  tests, 
the  spacing  between  adjacent  disks  was  5.08  cm  (two  disk  diameters),  as  illustrated  in  Fig¬ 
ure  2.  Based  on  the  information  gained  from  the  single  disk  simulations,  the  multiple-disk 
calculations  used  a  smoothing  length  of  0.03175  cm,  with  one  particle  per  smoothing  length 
yielding  ten  particles  across  the  thickness  of  each  disk. 
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Figure  6.  Predicted  penetration  channel  for  single-disk  impact,  particle  size  of: 
(a.)  0.0635  cm,  (b.)  0.03175  cm,  (c.)  0.02117  cm,  and  (d.)  0.015875  cm. 
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Figure  7.  Single-disk  penetration  depth  as  a  function  of  particle  size. 


Table  2.  Calculation  Characteristics. 


Particle  Size  (cm 

0.0635 

0.03175 

0.0217 

0.015875 

Particles  Across  Disk  Thickness 

5 

10 

15 

20 

Particles  in  Disk 

100 

400 

900 

1600 

Paxticles  in  Target 

24649 

99225 

222784 

396900 

Total  Particles 

24749 

99625 

223684 

398500 

Cycles^ 

2527 

5135 

7766 

10317 

Penetration  Depth  (cm) 

1.34 

1.43 

1.43 

1.43 

Computer  Time'*'  (hours) 

1.0 

8.2 

27.3 

66.5 

I  100  jUS  simulated  time 
t  195  MHz  Silicon  Graphics  RIOOOO  processor 
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Figure  8  illustrates  the  penetration  of  the  two  disks  into  the  target  at  simulated  times 
of  32  /is,  35  fis,  42  iJ,s,  and  100  fis.  At  32  fis,  the  initial  contact  between  the  two  disks 
occurs  at  the  outer  radius  of  the  second  disk.  Before  this  initial  contact  between  the  two 
disks,  the  calculation  was  identical  to  the  single-disk  case.  This  initial  contact  at  the  outer 
radius  initiates  a  set  of  radially  propagating  compression  waves  that  will  ultimately  converge 
at  the  centerline  of  the  disk.  By  35  fxs,  the  face  of  the  second  disk  at  the  centerline  has 
contacted  the  first  disk.  The  convergence  of  the  radially  propagating  compression  waves  at 
the  centerline  of  the  disk,  coupled  with  the  waves  generated  by  the  impact  of  the  face  of  the 
disk  at  the  center,  results  in  a  strong  expansion  wave  that  is  responsible  for  failure  of  the 
disk  material  near  the  center.  The  image  in  Figure  8  for  35  ixs  provides  the  initial  indication 
of  material  failure  as  the  particles  near  the  center  of  the  disk  are  beginning  to  diverge  at 
this  time,  sending  a  jet  of  failed  disk  material  in  the  opposite  direction  of  the  disk  motion. 
At  42  yus,  the  void  in  the  center  of  the  disk  has  grown  to  its  maximum  size  with  the  bulk  of 
the  disk  material  still  carrying  momentum  in  the  direction  of  the  target.  By  the  end  of  the 
simulation  at  100  yus,  this  momentum  has  caused  the  void  at  the  center  to  partially  close, 
with  a  small  void  remaining  at  the  center. 

To  further  illustrate  the  cause  of  failure  of  the  second  disk,  the  pressure  obtained  from 
the  tracer  at  the  center  of  the  leading  face  of  this  disk  is  presented  in  Figure  9.  This  figure 
has  two  history  plots.  The  upper  plot  shows  the  full  100-/l(s  range  of  tracer  data.  At  the 
time  of  impact  of  the  second  disk,  the  pressure  reaches  a  peak  of  approximately  4  Mbar. 
The  lower  plot  provides  a  magnified  view  of  the  data  centered  about  the  time  of  impact  of 
the  second  disk.  This  plot  shows  that  the  initial  pressure  peak  in  the  upper  plot  occurs  at 
approximately  33  yus  and  has  a  duration  of  approximately  1  yus.  After  this  positive  pressure 
phcLse,  the  pressure  becomes  negative  and  quickly  reaches  the  tensile  limit  of  the  tungsten 
alloy  at  -30  kbar.  Because  the  material  is  not  permitted  to  support  tension  beyond  this 
threshold  pressure,  the  curve  is  flat  at  -30  kbar  until  the  particle  again  sees  compression.  A 
second  negative  phase  occurs  between  35  yus  and  36  fis,  which  again  exceeds  the  tensile  limit 
of  the  disk  material. 

Of  the  set  of  two  disk  experiments  that  were  conducted,  only  one  had  sufficiently  low 
yaw  angle  of  both  disks  and  axial  alignment  of  disks  to  qualify  it  as  an  axisymmetric  event. 
For  this  experiment,  the  final  penetration  depth  into  the  target  material  was  2.35  cm.  The 
calculation  produced  a  total  penetration  depth  of  2.78  cm,  which  is  approximately  18% 
greater  than  the  experiment.  Since  there  is  only  one  experimental  data  point,  it  is  impossible 
to  determine  the  error  in  the  experimental  data  for  this  two-disk  configuration.  However, 
one  would  expect  the  experimental  scatter  to  be  at  least  as  large  as  that  for  the  single-disk 
case. 

The  second  disk  was  recovered  from  the  experiment  and  a  photograph  of  this  disk  is 
provided  in  Figure  10.  The  impact  face  of  the  disk  is  facing  the  camera  and  the  photograph 
shows  a  hole  in  the  disk  as  predicted  by  the  calculation.  One  can  see  that  the  hole  in  the 
recovered  disk  is  not  exactly  in  the  center  of  the  disk  and  that  the  disk  is  not  perfectly  round. 
These  characteristics  are  considered  to  be  a  result  of  small  asymmetries  in  the  experiment. 
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Figure  8.  Two-disk  impact  at  32  /us,  35  /is,  42  fis,  and  100  /us  (/i=0.03175  cm). 
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Figure  9.  Pressure  at  center  of  leading  face  of  second  disk. 
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Figure  10.  Recovered  projectile  from  two-disk  experiment. 

The  width  of  the  recovered  disk  varies  approximately  from  2.0  cm  (measured  vertically 
on  the  photograph)  to  2.3  cm  (measured  horizontally  on  the  photograph).  The  final  state 
of  the  calculation,  as  illustrated  in  the  100-^s  frame  of  Figure  8,  shows  that  the  second  disk 
has  some  eroded  material  deposited  along  the  wall  of  the  crater  and  the  eroded  material  is 
separated  from  the  main  body  of  the  disk.  The  main  body  has  the  void  region  at  the  center 
and  a  radius  of  approximately  2.4  cm,  which  compares  well  to  the  size  of  the  recovered  disk. 

The  final  axisymmetric  calculation  run  in  this  study  was  a  simulation  of  a  four-disk 
impact  experiment.  As  in  the  two-disk  case,  the  separation  distance  between  adjacent  disks 
was  5.08  cm,  and  the  smoothing  length  was  0.03175  cm,  with  one  particle  per  smoothing 
length.  Lagrangian  tracers  were  placed  on  the  impact  face  of  the  target  and  on  each  face 
of  the  disks  along  the  axis  of  symmetry.  The  calculation  was  run  to  a  simulated  time  of 
150  fis,  and  a  solution  image  was  saved  at  l-/xs  intervals  for  the  purpose  of  post-processing 
the  results. 

Figure  11  provides  material  plots  to  illustrate  four  of  the  most  significant  events  in  the 
simulation.  In  the  two-disk  case,  it  was  shown  that  the  second  disk  experienced  failure  at  the 
center  upon  impact  with  the  first  disk.  This  failure  sent  a  jet  of  disk  material  backwards,  in 
the  opposite  direction  of  the  disk  motion.  In  the  four-disk  simulation,  the  third  disk  initially 
encounters  this  debris  jet  of  material  after  45  /us  and  is  severely  damaged  by  the  debris  jet 
before  it  begins  to  interact  with  the  target  at  the  bottom  of  the  crater  created  by  the  first 
two  disks. 

The  material  plot  at  63  fis  shows  a  failed  region  in  the  center  of  the  third  disk  as  it 
makes  initial  contact  with  the  material  at  the  bottom  of  the  crater.  This  damage  decreases 
the  ability  of  the  third  disk  to  penetrate  the  target.  Some  of  the  material  that  caused  the 
failure  of  the  third  disk  at  the  centerline,  as  well  as  some  failed  material  from  the  third  disk 
itself,  continue  to  produce  a  debris  jet  that  encounters  the  fourth  disk  as  it  approaches  the 
bottom  of  the  crater.  The  material  plot  at  93  fis  shows  slight  damage  to  the  fourth  disk 
before  impact  on  the  crater.  The  large  void  region  near  the  centerline  left  by  the  failure  of 
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Figure  11.  Four-disk  impact  at  45  yus,  63  /j,s,  93  yus,  and  150  yus  (/i=0.03175  cm). 
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the  third  disk  causes  the  fourth  disk  to  fold  upward  near  the  outer  radius,  focusing  the  mass 
toward  the  centerline.  The  material  plot  at  150  fis  in  Figure  11  shows  the  final  state  of  the 
calculation  with  most  of  the  material  of  the  first  three  disks  deposited  along  the  walls  of 
the  crater  and  most  of  the  material  from  the  fourth  disk  at  the  bottom  of  the  crater  in  the 
center. 

The  position  history  of  the  Lagrangian  tracer  on  the  impact  face  of  the  target  can  be 
seen  in  Figure  12,  which  shows  four  distinct  impacts  by  the  individual  disks.  Immediately 
before  the  arrival  of  the  second  disk,  the  first  disk  had  penetrated  1.40  cm  into  the  target. 
When  the  third  disk  arrived,  the  penetration  of  the  second  disk  was  almost  complete  and 
had  contributed  an  additional  1.41  cm.  The  same  is  true  for  the  third  disk;  yet,  it  only 
contributes  0.79  cm  to  the  overall  penetration  depth.  Finally,  the  fourth  disk  adds  another 
1.36  cm,  to  bring  the  total  penetration  depth  to  4.96  cm.  Analysis  of  the  data  in  Figure  12 
shows  that  the  first,  second,  and  fourth  disks  provide  almost  equal  contribution  to  the  overall 
penetration,  while  contribution  from  the  third  disk  was  less  than  60%  of  any  of  the  other 
three.  The  performance  degradation  is  a  result  of  the  interaction  of  the  third  disk  with  the 
failed  material  (i.e.,  debris  jet)  of  the  second  disk  and  subsequent  centerline  failure. 


Figure  12.  Tracer  position  history  of  four-disk  impact  (A=0.03175  cm). 

Like  the  two-disk  case,  there  exists  only  one  experiment  using  four  disks  with  disk 
alignment  and  yaw  qualities  that  can  be  characterized  as  an  cixisymmetric  configuration. 
For  this  experiment,  the  total  penetration  depth  of  the  four  disks  into  the  semi-infinite  RHA 
target  was  4.45  cm,  approximately  11.5%  less  than  that  predicted  by  the  calculation.  As 
with  the  two-disk  experiment,  the  scatter  in  the  experimental  data  for  this  test  configuration 
is  not  known. 
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5.  THREE-DIMENSIONAL,  SINGLE  DISK  SIMULATIONS 


To  evaluate  the  capabilities  of  the  3D  implementation  of  the  SPH  code,  a  set  of  cal¬ 
culations  was  performed  to  simulate  the  single-disk  impact  tests.  As  in  the  case  of  the 
axisymmetric  calculations,  a  smoothing  length  resolution  study  was  performed  to  assess  the 
ability  of  the  3D  code  to  converge  on  a  final  penetration  depth  as  the  problem  resolution  was 
improved.  This  study  would  have  the  benefits  of  allowing  for  a  direct  comparison  between 
the  results  of  geometrically  identical  2D  and  3D  calculations.  It  would  also  reveal  any  addi¬ 
tional  modeling  challenges  that  exist  in  3D  that  may  not  be  present  in  2D  and  would  provide 
data  about  the  amount  of  computer  resources  (time,  memory,  disk  storage)  that  would  be 
required  to  perform  a  realistic  simulation. 

The  primary  difference  between  the  3D  models  used  in  this  set  of  calculations  and  those 
used  in  the  axisymmetric  cases  is  the  definition  of  particles  in  the  target.  The  axiymmetric 
models  that  were  used  to  simulate  the  single-disk  experiments  modeled  the  target  as  a  solid 
block  with  a  radius  of  10  cm  and  a  depth  of  10  cm,  with  a  set  of  uniformly  sized  particles 
distributed  throughout  the  target.  Similarly,  the  3D  model  represented  the  target  material 
as  a  solid  block  with  a  half  width  (distance  from  impact  point  to  outer  edge)  of  12.7  cm  and 
a  depth  of  12.7  cm  (ten  disk  radii).  However,  the  distribution  of  a  uniform  set  of  particles 
throughout  this  target  would  have  resulted  in  an  unacceptably  large  number  of  particles  in 
the  target,  resulting  in  very  large  computational  requirements. 

To  reduce  the  number  of  particles  in  the  target,  the  3D  computational  model  was  defined 
so  that  target  particles  identical  in  size  to  the  particles  in  the  penetrator  were  used  in  the 
region  where  the  interaction  between  the  two  bodies  was  to  occur.  Successive  layers  of 
increasing  particle  size  were  then  used  to  transition  from  this  central  region  of  interest  to 
the  boundaries  of  the  computation.  Figure  13  illustrates  the  particle  size  variation  in  the 
target  for  the  3D  model. 

In  the  figure,  the  target  particles  are  shaded  and  sized  according  to  the  value  of  the 
smoothing  length,  h.  Six  regions  of  constant  smoothing  length  are  used  to  vary  the  particle 
size  across  the  dimensions  of  the  target.  These  regions  are  identified  by  the  hi  through  he 
labels  in  Figure  13.  The  ratio  of  particle  sizes  in  two  adjacent  regions  is  1.75.  A  listing  of 
all  the  smoothing  lengths  used  in  the  3D  calculations  is  provided  in  Table  3. 


Table  3.  Smoothing  Lengths  used  in  3D  Target  Models. 


Particles 
Across  Disk 
Thickness 

Smoothing  Length  (cm) 

Total 
Particles 
in  Model 

hi 

h2 

hs 

h4 

hs 

he 

5 

0.06350 

0.11113 

0.19447 

0.34032 

0.59556 

1.04223 

187894 

6 

0.05292 

0.09260 

0.16206 

0.28360 

0.49630 

0.86853 

309554 

7 

0.04536 

0.07938 

0.13891 

0.24309 

0.42540 

0.74445 

488568 

8 

0.03969 

0.06945 

0.12154 

0.21270 

0.37223 

0.65139 

731090 

9 

0.03528 

0.06174 

0.10804 

0.18907 

0.33087 

0.57902 

1034672 

10 

0.03175 

0.05556 

0.09723 

0.17016 

0.29778 

0.52112 

1442442 
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Figure  13.  Target  particle  size  variation  for  3D  model. 

Calculations  using  the  models  described  in  Table  3  were  run  for  a  simulated  time  of 
40  /is.  This  simulation  time  is  based  on  the  previous  2D  single-disk  impact  simulations  that 
showed  the  penetration  event  to  be  complete  by  40  //s.  Because  the  code  employs  an  explicit 
solution  method,  the  time  required  to  complete  a  3D  simulation  is  proportional  to  the  fourth 
power  of  the  inverse  of  the  particle  size.  Thus,  reducing  the  particle  size  by  a  factor  of  two 
will  increase  the  run  time  by  a  factor  of  16.  The  same  is  true  for  explicit  grid-based  methods. 
In  the  case  of  the  3D  calcualtions  presented  here,  the  case  using  a  minimum  particle  size  of 
0.06350  cm  completed  in  11  hours.  However,  the  calculation  that  used  a  minimum  particle 
size  of  0.03175  cm  required  186  hours  to  complete.  This  is  consistent  with  the  fourth  power 
estimate  from  the  first  calculation.  The  run  times  of  the  other  calculations  in  the  set  were 
likewise  consistent  with  the  estimated  times. 

Figure  7  compares  the  penetration  depths  of  the  axisymmetric  calculations  to  the  exper¬ 
imental  data.  Figure  14  is  a  reproduction  of  that  figure,  with  the  addition  of  the  penetration 
depths  from  the  3D  calculations.  As  can  be  seen  in  the  figure,  the  3D  calculations  pro¬ 
duce  a  different  convergence  behavior  than  the  axisymmetric  code.  For  all  particle  sizes, 
the  final  penetration  depth  of  the  3D  code  is  less  than  its  axisymmetric  equivalent.  Like 
the  axisymmetric  calculations,  the  calculation  using  the  largest  particles  produced  the  least 
penetration  depth.  The  penetration  depth  increases  with  problem  resolution  up  to  a  particle 
size  of  0.03969  cm.  This  calculation  is  the  only  one  in  the  set  for  which  the  final  penetra¬ 
tion  depth  falls  within  the  range  of  uncertainty  of  the  experimental  data.  For  the  two  cases 
with  smaller  particles  than  this  case,  the  final  penetration  depths  were  slightly  lower.  Ad¬ 
ditional  calculations  with  even  smaller  particles  may  show  whether  or  not  true  convergence 
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was  obtained,  but  the  run  times  of  such  cases  would  be  prohibitively  long.  One  possible 
explanation  for  the  under-prediction  of  the  penetration  depth  of  the  3D  code,  as  compared 
to  the  axisymmetric  calculations,  may  be  the  variation  of  particle  sizes  in  the  target.  The 
interaction  of  particles  of  different  sizes  is  the  subject  of  the  next  section  of  this  report. 


Figure  14.  Single-disk  penetration  depth  as  a  function  of  particle  size  (2D  and  3D  results). 

6.  INTERACTION  OF  REGIONS  OF  DIFFERENT  PARTICLE  SIZES 


The  3D  calculations  described  in  the  previous  section  had  the  target  material  modeled 
with  multiple  sets  of  particles,  gradually  increasing  in  size  with  increasing  distance  from 
the  region  of  interest.  As  noted  earlier,  the  ability  of  the  3D  code  to  converge  on  a  final 
penetration  depth  with  increasing  problem  resolution  differed  significantly  from  the  results 
of  the  axisymmetric  simulations.  One  possible  explanation  for  this  difference  is  this  use  of 
varying  particle  sizes  in  the  target  model. 

To  further  study  the  effect  of  different  particle  sizes,  a  set  of  calculations  was  run  in  which 
the  size  of  the  particles  in  the  target  differed  from  the  size  of  the  particles  in  the  penetrator. 
For  all  these  single-disk  calculations,  the  penetrator  had  a  constant  smoothing  length  of 
0.03175  cm,  using  one  particle  per  smoothing  length.  For  each  calculation,  the  smoothing 
length  of  the  target  was  varied,  using  one  particle  per  smoothing  length.  Calculations  were 
run  in  which  the  relative  sizes  of  the  target  particles,  with  respect  to  the  penetrator  particles, 
were  0.50,  0.65,  0.80, 1.25,  1.30,  and  2.00.  These  cases  resulted  in  target  smoothing  lengths  of 
0.01588  cm,  0.02064  cm,  0.02540  cm,  0.03969  cm,  0.04128  cm,  and  0.06350  cm,  respectively. 
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The  case  in  which  the  target  particles  are  1.30  times  the  size  of  the  penetrator  is  special 
in  that  this  configuration  produces  target  and  penetrator  particles  of  nearly  identical  mass. 
The  calculations  were  run  to  40  yus,  and  the  solutions  were  compared  to  the  baseline  case  in 
which  the  target  and  penetrator  had  equal  sized  particles  (a  relative  target  particle  size  of 
1.00). 

Material  plots  at  40  ij,s  for  the  cases  using  particle  size  ratios  of  0.50,  1.00,  1.30,  and 
2.00  are  shown  in  Figure  15.  These  plots  show  the  dramatic  effect  of  the  variation  of  the 
particle  size  in  the  target  relative  to  that  in  the  penetrator.  The  material  plot  in  Figure  15a 
shows  the  result  for  the  case  in  which  the  penetrator  particles  are  twice  the  size  of  the  target 
particles.  For  this  case  and  the  others  in  which  the  ratio  of  target  to  penetrator  particle  size 
was  less  than  unity,  the  disk  penetrated  much  deeper  into  the  target  than  the  experiment.  In 
fact,  at  40  yus,  the  penetration  is  still  not  complete  for  these  cases.  Defining  the  penetrator 
particles  to  be  larger  than  those  in  the  target  tends  to  make  the  penetrator  more  rigid  and 
increases  its  ability  to  penetrate  the  target. 

The  material  plot  in  Figure  15b  shows  the  result  of  the  baseline  calculation  in  which  the 
penetrator  and  the  target  were  defined  with  identical  particle  sizes  (ratio  =  1.00).  Figure  15c 
and  15d  provide  the  material  plots  for  two  of  the  three  cases  in  which  the  target  particles 
were  larger  than  those  in  the  penetrator.  For  each  calculation  of  this  type  (ratio  >  1.00)  the 
target  material  behaved  in  a  more  rigid  manner  than  in  the  baseline  case.  The  most  extreme 
case  was  the  one  in  which  the  target  particles  were  twice  the  size  of  the  penetrator  particles. 
In  this  case,  the  penetrator  leaves  only  a  small  indentation  in  the  target  and  then  bounces 
off  the  target  and  breaks  into  pieces. 

These  results  help  to  explain  the  difference  found  between  the  convergence  properties 
of  the  2D  and  3D  solutions.  It  is  possible  that  if  the  3D  calculations  described  earlier  had 
been  run  with  a  constant  particle  size,  instead  of  varying  the  particle  size  in  the  target,  the 
3D  code  may  have  produced  greater  final  penetrations,  bringing  them  in  closer  agreement 
with  the  2D  results  and  the  experimental  data.  Further  examination  of  the  formulation  of 
the  SPH  algorithm  is  required  to  determine  the  cause  of  the  huge  variations  between  these 
calculations. 


7.  VARIATION  OF  SMOOTHING  LENGTH  RETAINING  CONSTANT 
PARTICLE  SIZE 

All  the  calculations  described  so  far  were  set  up  to  define  particle  sizes  as  being  equal 
to  the  smoothing  length.  Put  another  way,  each  calculation  used  only  one  particle  per 
smoothing  length.  It  is  not  necessary,  however,  to  set  up  problems  in  this  manner.  One 
possible  variation  on  problem  definition  would  be  to  increase  both  the  smoothing  length 
and  the  number  of  particles  per  smoothing  length,  thereby  retaining  the  particle  size  and 
total  number  of  particles  in  the  simulation.  Increasing  the  smoothing  length  would  have  two 
primary  effects  on  the  solution.  First,  the  number  of  neighbors  of  a  particular  particle  would 
increase,  thereby  increasing  the  number  of  particles  used  in  updating  the  state  variables  at 
every  time  step.  This  would  create  more  work  to  be  done  at  every  time  step  for  each  particle. 
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Figure  15.  Single-disk  impact  at  40  /xs  for  particle  size  ratios  of:  (a.)  0.50,  (b.)  1.00,  (c.)  1.30, 
and  (d.)  2.00. 
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For  2D  calculations,  the  number  of  neighboring  particles  used  in  updating  the  state  variables 
is  proportional  to  the  square  of  the  smoothing  length,  resulting  in  a  workload  increase  that 
is  proportional  to  h'^.  Second,  since  the  size  of  the  time  step  is  directly  proportional  to  the 
smoothing  length,  the  time  step  would  increase,  resulting  in  fewer  iteration  cycles  required 
to  complete  the  simulation,  resulting  in  a  decrease  in  workload  that  is  proportional  to  h. 
Therefore,  for  a  set  of  2D  calculations  of  this  type,  the  time  required  to  complete  a  simulation 
should  be  proportional  to  the  smoothing  length,  h.  Using  this  logic,  it  can  be  seen  that 
the  time  required  for  a  set  of  3D  calculations  would  be  proportional  to  the  square  of  the 
smoothing  length. 

A  final  set  of  axisymmetric  calculations  was  performed  by  taking  this  approach  to  model 
the  single-disk  impact  experiments.  A  constant  particle  size  of  0.03175  was  employed,  re¬ 
sulting  in  a  total  of  99,625  particles  in  each  calculation.  Each  calculation  was  run  to  a 
simulated  time  of  40  fis.  As  in  the  other  calculations,  Lagrangian  tracer  particles  were  used 
to  gather  position  histories  of  the  penetrator  and  target  materials.  The  characteristics  of 
these  calculations  and  their  resulting  final  penetration  depths  are  provided  in  Table  4. 


Table  4.  Calculations  in  Smoothing  Length  Variation  Study. 


h 

(cm) 

Number  of 
Particles 
per  h 

Iteration 

Cycles 

Computer 

Time 

(hours) 

Penetration 

Depth 

(cm) 

0.03175 

1 

0.06350 

2 

HUH 

0.12700 

4 

581 

11.0 

1.21 

0.25400 

8 

286 

22.0 

1.08 

0.50800 

16 

160 

50.7 

0.60 

Table  4  shows  that  the  estimated  workload  associated  with  these  calculations  closely 
follows  the  expectation;  the  number  of  iteration  cycles  decreased  in  proportion  to  the  in¬ 
crease  in  smoothing  length,  and  the  total  computer  time  increased  proportionally  with  the 
smoothing  length.  The  most  dramatic  result  of  this  set  of  calculations  was  the  effect  of 
smoothing  length  on  the  resulting  simulation.  Table  4  shows  that  as  the  smoothing  length 
was  increased,  the  penetration  depth  decreased  significantly.  This  decrease  in  penetration 
depth  initially  suggests  that  the  increase  in  smoothing  length  tends  to  artificially  increase 
the  rigidity  of  the  materials  and  therefore  inhibit  penetration.  Further  examination  of  the 
results  shows  that  the  calculations  differ  not  only  in  final  penetration  depth  but  in  the  basic 
characteristics  of  the  impact  event.  Figure  16  shows  the  y-position  (axial)  histories  of  the 
tracer  particle  located  at  the  center  of  the  target  for  each  of  the  calculations.  The  solid  line 
presents  the  results  of  the  baseline  calculation  using  one  particle  per  smoothing  length.  A 
detailed  description  of  the  penetration  event  for  this  calculation  was  provided  earlier  in  this 
report.  The  figure  shows  that  calculations  employing  more  than  one  particle  per  smoothing 
length  produced  entirely  different  penetration  histories  from  the  baseline  calculation.  The 
final  penetration  depth  produced  by  the  baseline  calculation  is  in  closer  agreement  with 
the  experimental  data  than  those  calculations  using  more  than  one  particle  per  smoothing 
length.  Furthermore,  calculations  with  other  grid-based  codes  produced  results  simular  to 
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those  created  by  the  baseline  calculation.  These  findings  suggest  that  it  is  best  to  run  Magi 
calculations  of  high  velocity  impact  with  one  particle  per  smoothing  length.  Further  study 
of  the  numerical  formulation  of  the  SPH  method  is  necessary  to  determine  the  cause  of  the 
significant  differences  between  these  calculations. 


Figure  16.  Penetration  Histories  from  Smoothing  Length  Variation  Study. 


8.  SUMMARY  AND  RECOMMENDATIONS 


This  report  has  documented  the  simulation  of  a  set  of  disk-shaped  penetrator  impact 
experiments  with  the  SPH  code  Magi.  The  2D  implementation  of  the  code  was  tested 
extensively  to  reveal  some  of  the  strengths  and  weaknesses  of  the  SPH  method  for  simulating 
high  velocity  impact  problems.  Particle  resolution  studies  were  performed  to  determine  the 
ability  of  the  2D  and  3D  implementations  of  the  code  to  converge  on  a  final  penetration 
depth.  These  calculations  revealed  that  approximately  10  particles  are  required  across  any 
critical  dimension  to  properly  capture  the  transient  compression  and  expansion  waves  that 
occur  during  such  impact  events.  In  the  case  of  the  disk  impact  cases  examined  here, 
this  critical  dimension  is  the  thickness  of  the  penetrator.  As  the  particle  size  was  refined, 
the  calculational  results  fell  within  the  range  of  uncertainty  of  the  experimental  data  for 
a  single-disk  impact  event,  6.5%  below  the  mean  of  the  experimental  mean  penetration 
depth.  Calculations  were  also  performed  for  cases  involving  two  and  four  axially  aligned 
disk  projectiles  with  penetration  depths  predicted  to  within  11.5%  of  experimental  results. 

Additional  sets  of  calculations  were  performed  to  determine  the  effects  of  increasing  the 
number  of  particles  per  smoothing  length  and  the  interaction  of  different  sized  particles. 
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These  calculations  produced  results  that  differed  drastically  from  the  experimental  data. 
Overall,  the  best  results  were  obtained  by  using  a  contant  particle  size  throughout  the 
computational  domain  with  one  particle  per  smoothing  length.  The  SPH  method  and  specific 
Magi  programming  must  be  studied  in  greater  detail  to  understand  the  possible  causes  of 
these  problems. 

Of  course,  these  disk  impact  simulations  alone  are  not  sufficient  to  fully  evaluate  the 
viability  of  SPH  in  armor/anti-armor  applications.  However,  they  do  provide  a  positive  ini¬ 
tial  step  in  assessing  the  dominant  computational  characteristics  of  this  method.  Additional 
calculations,  particulary  in  3D,  over  a  wide  variety  of  configurations  of  interest  are  required 
to  complete  the  evaluation.  Experimental  data  exist  for  yawed,  oblique  long-rod  penetrator 
impact,  and  the  detonation/formation  of  an  explosively  formed  projectile.  Further  investi¬ 
gation  of  SPH  through  its  application  to  problems  of  this  type  is  essential  to  formulating  an 
accurate  assessment  of  the  technology. 
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